Load Library
library(spdep)
library(tidyverse)
library(readr)
library(readxl)
library(car)
library(sf)
library(stplanr)
library(terra)
library(tmap)
library(leaflet)
library(tseries)
library(dplyr)
library(sfweight)
library(ggplot2)
library(rgdal)
Load data
dataset <- read_excel("E:/FOLDER KULIAH/STIS/SEMESTER 7/SKRIPSI/dataset_skripsi_revisi.xlsx")
# Check missing data
sapply(dataset, function(x) sum(is.na(dataset)))
## No KodeKecamatan NamaKecamatan
## 0 0 0
## PEND_LGSG ntl_median lst_median
## 0 0 0
## bui_median co_median pm2.5_median
## 0 0 0
## elevation_median aridity_index_median rwi_weighted
## 0 0 0
## jlh_klg_listrik_pln jlh_fas_pddk jlh_fas_ksht
## 0 0 0
## jlh_imk jlh_sentra_industri jlh_sarana_eko
## 0 0 0
## jlh_tmpt_ibdh kepadatan_penduduk
## 0 0
str(dataset)
## tibble [572 x 20] (S3: tbl_df/tbl/data.frame)
## $ No : num [1:572] 1 2 3 4 5 6 7 8 9 10 ...
## $ KodeKecamatan : num [1:572] 3301010 3301020 3301030 3301040 3301050 ...
## $ NamaKecamatan : chr [1:572] "Dayeuhluhur" "Wanareja" "Majenang" "Cimanggu" ...
## $ PEND_LGSG : num [1:572] 23.9 16.1 19.8 20.6 25.5 ...
## $ ntl_median : num [1:572] 2.06 2.06 1.97 2.55 2.75 ...
## $ lst_median : num [1:572] 27.3 28.6 27.1 28.6 28.8 ...
## $ bui_median : num [1:572] -1.09 -1.04 -1.04 -1.02 -1.07 ...
## $ co_median : num [1:572] 0.0296 0.0307 0.0301 0.032 0.0293 ...
## $ pm2.5_median : num [1:572] 40.4 35.1 35.7 33.1 30.3 ...
## $ elevation_median : num [1:572] 347 113 214 155 112 52 15 6 3 14 ...
## $ aridity_index_median: num [1:572] 2.13 1.8 2.04 2.03 1.88 ...
## $ rwi_weighted : num [1:572] -0.066 0.4655 0.4854 0.2715 0.0747 ...
## $ jlh_klg_listrik_pln : num [1:572] 20142 40280 50842 38141 28380 ...
## $ jlh_fas_pddk : num [1:572] 119 191 277 180 166 149 136 194 85 208 ...
## $ jlh_fas_ksht : num [1:572] 40 75 82 65 44 42 53 48 32 68 ...
## $ jlh_imk : num [1:572] 157 357 490 675 162 ...
## $ jlh_sentra_industri : num [1:572] 0 0 2 2 4 0 0 0 0 0 ...
## $ jlh_sarana_eko : num [1:572] 513 839 934 521 503 ...
## $ jlh_tmpt_ibdh : num [1:572] 356 585 684 573 492 375 280 421 266 539 ...
## $ kepadatan_penduduk : num [1:572] 256 542 842 632 639 ...
# Filter variabel
data_x <- dataset[, -c(1,2,3,4)]
data_y <- dataset[, 4]
# Analisis Deskriptif
summary(data_x)
## ntl_median lst_median bui_median co_median
## Min. : 1.373 Min. :22.64 Min. :-1.20879 Min. :0.01706
## 1st Qu.: 2.216 1st Qu.:28.64 1st Qu.:-0.98658 1st Qu.:0.02682
## Median : 2.639 Median :30.57 Median :-0.75872 Median :0.02866
## Mean : 3.128 Mean :30.30 Mean :-0.77524 Mean :0.02882
## 3rd Qu.: 3.182 3rd Qu.:31.93 3rd Qu.:-0.62883 3rd Qu.:0.03111
## Max. :19.035 Max. :37.88 Max. : 0.03596 Max. :0.03821
## pm2.5_median elevation_median aridity_index_median rwi_weighted
## Min. :17.50 Min. : 1.0 Min. :0.8404 Min. :-0.4047
## 1st Qu.:32.10 1st Qu.: 25.0 1st Qu.:1.2586 1st Qu.: 0.2410
## Median :36.33 Median : 114.5 Median :1.4555 Median : 0.4299
## Mean :35.75 Mean : 246.5 Mean :1.5689 Mean : 0.4561
## 3rd Qu.:40.20 3rd Qu.: 331.0 3rd Qu.:1.8331 3rd Qu.: 0.6453
## Max. :49.40 Max. :1816.0 Max. :2.7325 Max. : 1.7157
## jlh_klg_listrik_pln jlh_fas_pddk jlh_fas_ksht jlh_imk
## Min. : 3516 Min. : 31.0 Min. : 11.00 Min. : 60.0
## 1st Qu.:15056 1st Qu.:104.0 1st Qu.: 38.00 1st Qu.: 301.8
## Median :19732 Median :133.5 Median : 50.00 Median : 558.5
## Mean :21474 Mean :142.6 Mean : 56.86 Mean : 883.8
## 3rd Qu.:25679 3rd Qu.:171.5 3rd Qu.: 69.00 3rd Qu.:1081.5
## Max. :61913 Max. :405.0 Max. :219.00 Max. :9152.0
## jlh_sentra_industri jlh_sarana_eko jlh_tmpt_ibdh kepadatan_penduduk
## Min. : 0.00 Min. : 28.0 Min. : 50.0 Min. : 117.0
## 1st Qu.: 0.00 1st Qu.: 556.2 1st Qu.:206.0 1st Qu.: 721.8
## Median : 1.00 Median : 789.5 Median :273.5 Median : 1121.5
## Mean : 2.57 Mean : 854.0 Mean :294.6 Mean : 1747.4
## 3rd Qu.: 3.00 3rd Qu.:1030.0 3rd Qu.:360.5 3rd Qu.: 1736.8
## Max. :54.00 Max. :3322.0 Max. :826.0 Max. :16094.0
summary(data_y)
## PEND_LGSG
## Min. : 8.309
## 1st Qu.:20.531
## Median :26.146
## Mean :27.604
## 3rd Qu.:32.614
## Max. :71.308
Uji Moran’s I
# moran.test
W <- read.gal("E:/GEOJSON/QC_MAT.gal", override.id=TRUE)
w <- nb2listw(W, glist=NULL, style="W", zero.policy=NULL)
# write.nb.gal(w$neighbours, "E:/GEOJSON/queen_cont.gal")
Y <- dataset$PEND_LGSG
result_moran <- moran.test(Y, w, randomisation=FALSE, zero.policy=TRUE, alternative="two.sided",
rank = FALSE, na.action=na.fail, spChk=NULL, adjust.n=TRUE)
result_moran
##
## Moran I test under normality
##
## data: Y
## weights: w
##
## Moran I statistic standard deviate = 7.084, p-value = 1.4e-12
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.1804663118 -0.0017513135 0.0006616349
Simulasi Monte Carlo
set.seed(123)
MC<- moran.mc(Y, w, nsim=599)
# View results (including p-value)
MC
##
## Monte-Carlo simulation of Moran I
##
## data: Y
## weights: w
## number of simulations + 1: 600
##
## statistic = 0.18047, observed rank = 600, p-value = 0.001667
## alternative hypothesis: greater
Local Morans’I
oid <- order(Y)
localMI <- localmoran(Y, w)
signifikansi_y <- rep(0,572)
localMI <- cbind(localMI, signifikansi_y)
localMI[,5] <- round(localMI[,5], 4)
for(i in 1:nrow(localMI)){
if(localMI[i,5] <= 0.05){
localMI[i,6] <- 1 # SIGNIFIKAN = 1
} else{
localMI[i,6] <- 0 # TIDAK SIGNIFIKAN = 0
}
}
# Jumlah kecamatan yg signifikan
sum(localMI[,6] == 1)
## [1] 73
# Jumlah kecamatan yg tidak signifikan
sum(localMI[,6] != 1)
## [1] 499
View(localMI)
Mapping Local Morans’I P-values
jateng.map <- st_read("E:/GEOJSON/join_new_data.shp",quiet=TRUE)
# MAPPING LOCAL MORAN'S I
listrik <- jateng.map$LISTRIK
listrik_localMI <- cbind(listrik, localMI)
# install.packages("tmap")
library(tmap)
# tm_shape(localMI1) +
# tm_fill(col = "Ii",
# style = "pretty",
# palette = "RdBu",
# title = "local moran statistics") +
# tm_borders(alpha = 0.5)
MoranLocalMap = cbind(listrik_localMI, jateng.map)
MoranLocalMap = st_as_sf(MoranLocalMap, crs=st_crs(MoranLocalMap$geometry))
current.mode <- tmap_mode("view") #switch to interactive map
## tmap mode set to interactive viewing
tm_shape(MoranLocalMap) +
tm_fill(col = "Pr.z....E.Ii..",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette = "-Blues",
title = "local Moran's I p-values")+
tm_borders(alpha = 0.5) +
tm_layout(title="Mapping Local Moran's I P-Value")
Mapping Local Morans’I values
tm_shape(MoranLocalMap) +
tm_fill(col = "Ii",
style = "pretty",
palette = "RdBu",
title = "local moran statistics (Ii)")+
tm_borders(alpha = 0.5) +
tm_layout(title="Mapping Local Moran's I Values")
## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
LISA Cluster MAP 4 quadrant
quadrant <- vector(mode = "numeric", length = nrow(localMI))
DV <- jateng.map$LISTRIK - mean(jateng.map$LISTRIK)
C_mI <- localMI[,1] - mean(localMI[,1])
signif <- 0.05
quadrant[DV >0 & C_mI>0] <- 4
quadrant[DV <0 & C_mI<0] <- 1
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI[,5]>signif] <- 0
MoranLocalMap$quadrant <- quadrant
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
# knitr::kable(head(MoranLocalMap, n=10))
tm_shape(MoranLocalMap) +
tm_fill(col = "quadrant", style = "cat", palette = colors[c(sort(unique(quadrant)))+1],
labels = clusters[c(sort(unique(quadrant)))+1], popup.vars = c("nmkec")) +
tm_view() +
tm_borders(alpha=0.5) +
tm_layout(title="LISA Map Cluster")